A FAMILY OF MIXED FINITE ELEMENT PAIRS WITH OPTIMAL 
GEOSTROPHIC BALANCE 

COLIN COTTER* 

Abstract. Wc introduce a family of mixed finite element pairs for use on geodesic grids and 
with adaptive mesh refinement for numerical weather prediction and ocean modelling. We prove 
that when these finite element pairs are applied to the linear rotating shallow water equations, the 
gcostrophically balanced states are exactly steady, which means that the numerical schemes do not 
introduce any spurious inertia-gravity waves; this makes these finite element pairs in some sense 
optimal for numerical weather prediction and ocean modelling applications. We further prove that 
these finite element pairs satisfy an inf-sup condition which means that they are free of spurious 
pressure modes which would pollute the numerical solution over the timescales required for large- 
scale geophysical applications. We then discuss the extension to incompressible Euler-Boussinesq 
equations with rotation, and show that for the linearised equations the balanced states are again 
exactly steady on arbitrary unstructured meshes. We also show that the discrete pressure Poisson 
equation resulting from these discretisations satisfies an optimal stencil property. All these properties 
make the discretisations in this family excellent candidates for numerical weather prediction and 
large-scale ocean modelling applications when unstructured grids are required. 
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1. Introduction. The aim of this paper is to introduce a new family of finite 
element pairs for discretising large-scale ocean and atmosphere flows, state theorems 
about their wave propagation properties, and to explain why these properties are 
important for numerical weather prediction and ocean modelling. Section [l . 1 1 provides 
some motivational background, section [Q| describes the balance properties of model 
equations that we aim to preserve exactly in the numerical discretisations, and section 
1.3| explains how these properties can be modified by numerical discretisations. Section 
2] introduces the family of finite element pairs, and the optimal balance property is 
proved in section [3] It is also shown, by proving an inf-sup condition, that the finite 
element pairs are free from spurious pressure modes. In section [4] these properties are 
extended to the case of three-dimensional rotating stratified incompressible flow which 
applies to non-hydrostatic ocean modelling. Finally, section [5] provides a summary 
and outlook to future developments. 

1.1. Background. An operational weather forecasting system combines obser- 
vational data (such as satellite images, and pressure measurements on the Earth's 
surface, for example) with a numerical model of the atmosphere (solved on state- 
of-the-art parallel computers) to produce weather forecasts. The numerical model 
consists of a mathematically-consistent discretisation of the equations of motion for 
the atmosphere (known as the dynamic core) together with schemes (known as param- 
eterisations) for representing subgrid scale physics such as convection and turbulence, 
and physical processes such as cloud formation and atmospheric chemistry. Most cur- 
rent state-of-the art numerical weather prediction (NWP) models, such as the MET 
Office Unified Model [7] , make use of a latitude- longitude grid to construct the numer- 
ical approximations to the fluid dynamical equations that form the dynamical core of 
the model. However, recently there has been growing interest in more general hori- 
zontal discretisation schemes constructed using triangles or hexagons. This is for two 
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reasons. Firstly, geodesic grids (which are obtained by iterative refinement of an icosa- 
hedron using triangles or the dual grid which comprises a combination of hexagons 
plus exactly 12 pentagons) provide a much more uniform coverage over the sphere, 
which has possible advantages for accurate representation of wave propagation and 
avoids very fine grid cells near to the North and South poles. A number of groups are 
now using geodesic grids for weather and climate models [20l[14j[22]. Secondly, trian- 
gles facilitate adaptive mesh refinement much more easily, allowing regional models to 
be nested seamlessly in a global model, and even allowing dynamic mesh refinement 
in which the mesh is dynamically adapted during a forecast. However, the introduc- 
tion of adaptively-refined triangular grids calls for the careful design of new numerical 
schemes which correctly represent the large scale geophysical balances so that model 
forecasts can be run over sufficiently long times (the current forecast window is about 
a week). In this paper we introduce a new family of numerical discretisation meth- 
ods on triangular grids which shall be shown to optimally represent these geophysical 
balances on arbitrary unstructured grids. We also extend these properties to three 
dimensional rotating stratified incompressible flow so that they may be applied to 
non-hydrostatic ocean modelling. 

1.2. Model equations and geostrophic states. As a model problem, we 
consider the shallow-water equations on an /-plane 

u t + (u ■ V)u + fu 1 - + g\7D = 0, u=(u x ,u 2 ), u 1 - = (-u 2 , Mi), (1.1) 



where u is the horizontal velocity, D = D + rj is the total layer depth, 77 is the 
perturbation layer thickness, D is the mean layer thickness, k is the unit vector in 
the z-direction, / is the (constant) Coriolis parameter and g is the acceleration due 
to gravity. The boundary conditions are 



where dfl denotes the boundary of the two-dimensional domain 17, and n is the normal 
to dQ. These equations model the dynamics of a layer of hydrostatic incompressible 
fluid with constant density with a free surface and columnar motion so that the 
horizontal velocity is independent of depth. This can be thought of as a simple model 
for the ocean, or for a single layer in the atmosphere. For a derivation of these 
equations, see |5T], for example. Numerical methods for NWP and ocean modelling 
are often developed in two dimensions using the rotating shallow-water equations. The 
methods can then be extended to the primitive equations (three-dimensional equations 
of rotating stratified hydrostatic flow, see [21]) by building a mesh in layers in which 
the two-dimensional discretisation is applied in the horizontal plane. Hence, the 
shallow- water equations provide a useful testbed for numerical schemes and a stepping 
stone to the primitive equations which already poses many of the key challenges for 
designing good numerical methods for NWP and ocean modelling. 

A key dimensionless number in geophysical fluid dynamics is the Rossby number, 
defined by 



D t + V ■ (uD) = 0, 



(1.2) 



u ■ n — on d£l, 



(1.3) 



Ro = U/fL, 



where U is a typical velocity scale, and L is a typical horizontal length scale; the 
Rossby number measures the relative importance of the acceleration and Coriolis 
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terms. Geophysical flow problems in the small Rossby number limit are concerned 
with the states which satisfy 

fu 1 + gVD fa 

so that the Coriolis term approximately balances the pressure gradient. This is the 
state of geostrophic balance. It has long been observed (see |15II16) . for example) that 
if the initial conditions for the state variables Uq(x) and Dq(x) are initialised near 
to geostrophic balance, then this state will be approximately preserved for very long 
times. This has observed to be the case both in the low Rossby number limit, and also 
in the 0(1) limit and various mechanisms have been proposed for the maintenance and 
breakdown of this balance in various parameter regimes [SJ [T7] ) . Large scale flows 
in the atmosphere and ocean (such as those associated with the global circulation that 
determines global weather and climate) are observed to be in a state of geostrophic 
balance and hence it is important that numerical schemes used for weather forecasting, 
ocean modelling and climate prediction can correctly represent these balances. 

The ability of a numerical scheme to represent geostrophic balance can be exam- 
ined by studying wave propagation properties. Choosing a flat topography so that D 
is a constant, the linearised shallow- water equations are 

ut + fu 1 - + gVT? = 0, (1.4) 
t] t + DV -u = Q. (1.5) 

If the velocity u and free surface elevation rj are chosen in a state of perfect geostrophic: 
balance so that 

fu 1 - + gVr) = 0, 
then the velocity divergence satisfies 

V • U = • = -T)y X + T) X y = 0, V 1 - = (-8y,d X ) 

and so we have a steady state 

u t = 0, rjt = 0. 

Hence, all gcostrophically-balanced states are steady for the linearised equations. In 
the nonlinear case, the solutions variables evolve through the nonlinear advection 
terms which are small in the low Rossby number limit, and hence the geostrophic 
balance is maintained as a form of adiabatic invariance due to the small parameter 
Ro. 

If periodic boundary conditions are chosen, then one can perform a dispersion 
analysis for solutions of equations ( 1.4|1.5 1 by substituting the ansatz 



resulting in the dispersion relation 

u; (iv 2 - f - gD (k 2 + I 2 )) = 0. (1.6) 



This equation has three roots, with the u) = root corresponding to the steady 
geostrophic state; this branch becomes the Rossby wave branch when the model is 
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extended to the (3-plane model in which the Coriolis parameter / varies in the y- 
direction so that / = / + j3y for constants /o and j3. The other two roots give rise 
to dispersive waves, known as inertia- gravity waves; the k = I = case is known as 
an inertial oscillation in which the fluid undergoes solid body motion with a flat free 
surface. Since all of the roots are real, the state of geostrophic balance is stable under 
small perturbations. 

1.3. Geostrophic states for numerical methods. It is crucial that numerical 
methods for equations ( 1.1|1.2 1 do not generate spurious inertia-gravity waves when 



the solution is near to geostrophic balance. This can typically occur if the wave- 
propagation properties of the numerical method (i.e. the numerical discretisation of 
the linearised equations) do not correctly represent this balance. Given a discretisation 
of the linear system ( 1.4|1.5 1 it is simple to check the evolution of geostrophically- 



balanced states under this discretisation. One constructs initial conditions which 
satisfy the discrete form of geostrophic balance, steps the variables forward in time, 
and inspects the variables to check that the steady state is approximately preserved. 
This analysis has been performed for various element pairs in [13] . In general, the 
discrete divergence of the velocity field will not be exactly zero, and the remainder 
due to numerical discretisation errors will lead to oscillations. Whether the numerical 
method is suitable for computing the evolution of gcostrophically-balanced states 
over long time intervals {i.e. suitable for weather forecasting or ocean modelling) 
depends on how these errors behave. If one computes the numerical dispersion relation 



(i.e. the numerical analogue of equation (1.6| for the method (see for example the 



calculations in (23l [12]) then it is possible to divide the eigenmodes of the system into 
geostrophic modes which converge to the geostrophic states as the mesh edge-lengths 
converge to zero, and inertia-gravity modes which converge to the inertia-gravity 
waves. If the numerical discretisation errors in the divergence of the balanced states 
are large and project onto the inertia-gravity modes, then large unbalanced dynamics 
will be apparent after a long time integration interval (such as the time interval 
that is relevant to weather forecasting). When solving the nonlinear equations these 
numerical errors are constantly generated by the nonlinear terms, resulting in the 
geostrophic component of the solution being polluted by spurious inertia-gravity waves 
which render the numerical scheme useless for NWP and global ocean modelling. 

In section [2] we present a family of mixed finite element pairs which have the opti- 
mal property that the geostrophically balanced states are exactly steady; this means 
that these numerical discretisations are in some sense optimal for geophysical fluid 
dynamics problems. This property is independent of the choice of mesh, which can be 
taken to be completely unstructured. This means that the numerical discretisations 
can be used to solve geophysical flow problems in the presence of mesh refinement 
and adaptivity. This is proved in theorem |3.3| in section [3] We shall also show that 
the finite element pairs satisfy an inf-sup condition which means that they are free 
from spurious pressure modes: eigenmodes which have very small discrete gradients 
of free surface elevation despite having a free surface which is not flat. The absence of 
these modes is also crucial for geophysical applications since they can be coupled to 
the physical modes through the nonlinear terms and eventually becoming as large as 
the physical solution; the existence of these modes prohibits the use of such numerical 
methods in weather forecasting and ocean modelling. This result is proved in theorem 
3.7 also in section [3| The proofs of these results are very simple and elegant, due to 



the geometric embedding conditions that define the family of discretisations. 
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2. Family of finite element pairs. In this section we introduce our family 
of finite element pairs, first by developing the general finite element formulation in 
section |2~T| and then by stating conditions which define our particular family in section 



2.1. Finite element formulation. In this subsection we develop the finite el- 
ement approximation to the linearised shallow-water equations by writing down the 
weak form of the equations and restricting the function spaces to chosen spaces of 
piecewise polynomials on a finite clement mesh. We start with the linearised shallow- 
water equations on an /-plane given in equations ( 1.4|1.5 I with boundary conditions 
given by equation ( 1.3 1. To obtain the weak form of the equations we multiply equa- 
tion (1.4) by a test function w and equation (1.5) by a test function (j) and integrate 



over the domain ft to obtain 
d 
df 



w 



udV + f j w 

o 

_d_ 

df 



u^dV = -g / w ■ V?? d V, 
Jn 

<j>T]dV = -D <f>V ■ udV. 



(2.1) 
(2.2) 



We then integrate equation (2.2 1 by parts, and make use of the boundary conditions 



(1.3 1 to obtain 



— [ w ■ udV + f [ w ■ u^dV = -g [ w-VrjdV, 
dt Jq J n Jq 



m 
_d_ 

dt. 



<prjdV = D / Vcf)-udV, 



(2.3) 
(2.4) 



which must hold for all test functions w £ i/ 1 (SI) and <fr £ L 2 (£l). This is the weak 
form of equations ( 1.4|1.5 1 which we shall discretise using the Galerkin finite element 
method. 

The Galerkin projection of equations ( 2.3|2.4 | is constructed by defining finite 
dimensional spaces for the numerical solution variables u s and r/ 6 , and the test func- 
tions w s and 4> s . We shall use a mixed finite element method (see [3] for an excellent 
general survey) , which means that one type of finite clement space shall be used for 
u s and w s , and a different type of finite element space shall be used for r] S and <p s . 

We shall begin by defining the possible finite element spaces in general terms, 
before going on to state conditions which define our particular family of discretisations. 

Definition 2.1 (Finite element mesh). Let the mesh M. be a set of non- 
overlapping polygons (elements) which completely cover the computational domain 
fl which has a elementwise polygonal boundary dfl. 

Definition 2.2 (Pressure space). Let H be a space of elementwise polynomials 
on J\4, of type and continuity to be specified. This is a general definition, but we 
note that we shall require at least C° conti nuity ac ross element boundaries, since we 
apply gradients to <j) S and r] S in equations ( 2.5|2.6 1. 

Definition 2.3 (Velocity space). Let V be a space of vectors of elementwise 
polynomials on M, of type and continuity to be specified (possibly differently to H). 
We note that we do not require any co ntinuity conditions for V, since gradients are 
not applied to u s and w s in equations ( 2.5^2.6 1 . 

Having defined H and V, we may now write down the Galerkin finite element 
method for equations ( 1.4||1.5 1, which is obtained by restricting the solution variables 
and the test functions to these finite dimensional spaces. 
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Definition 2.4 (Galerkin finite element method). u s (x,t) and rj S (x,t) are the 



semi-discrete solutions of the Galerkin finite element discretisation of (1.4 1.5) if 
u\-,t)&V, v 5 (;t)£H, Vie[0,T], 



and 



w s ■ u s d V + f [ w s ■ u 5± d V = -g [ w s -VrfdV, (2.5) 
Jn Jn 

/ 4> B T) s dV = D [ Vc/) S -u s dV, (2.6) 



for all test functions w s £V,(p s £H. This equations may be solved on a computer 
by expanding w s and u s in a basis for V, and tfi s and h s in a basis for H , which 
produces a matrix equation for the basis coefficients of u s and h . This equation may 
then be discretised in time using a suitable time integration method. 

2.2. Choice of finite element spaces. In defining the problem, it remains to 
select a particular choice spaces (V, H) (known as a finite element pair). In this paper, 
we discuss a large family of possible choices defined by the following condition: 

Definition 2.5 (Embedding conditions) . 

1. The operator V defined by the pointwise gradient 

q S (x) = Vh 5 {x) 

maps from H into V . 

2. The skew operator _L defined by the pointwise formula 

q s (x) = (u^x)) 1 - 

maps from V into itself. 

These conditions are most definitely not satisfied by all possible pairs (V,H), as 
illustrated by the following examples. 

Example 2.6 (Pl-Pl). The finite element pair known as Pl-Pl (which may be 
used for the shallow-water equations but requires stabilisation as described in I25f ) is 
defined as follows: 

• The mesh M. is composed of triangular elements. 

• H is the space of elementwise-linear functions h which are continuous across 
element boundaries. 

• V is the space of vector fields u with both of the Cartesian components 
(u s ,v s ) in H. 

Condition 1 of Definition \2. 5| is not satisfied by the Pl-Pl pair since gradients of func- 
tions in H are discontinuous across element boundaries. Condition 2 is satisfied since 
the same continuity conditions are required for normal and tangential components. 

Example 2.7 (RTO). The lowest order Raviart- Thomas I19\j velocity space 
(known as RTO) is constructed on a mesh M composed of triangular elements. It 
consists of elementwise constant vector fields which are constrained to have continu- 
ous normal components across element boundaries. RTO does not satisfy condition 2 
of Definition \2.5\ since the _L operator transforms vector fields with discontinuities in 
the tangential component (which are permitted in RTO) into vector fields with discon- 
tinuities in the normal component (which are not). 



Mixed elements with optimal geostrophic balance 



7 



We now describe some examples of finite clement pairs which do satisfy the con- 
ditions in Definition 12.51 

Example 2.8 (P0-P1). The finite element pair known as P0-P1 (applied to ocean 
modelling in \21$ , for example) is defined as follows: 

• The mesh M. is composed of triangular elements. 

• H is the space of elementwise-linear functions h s which are continuous across 
element boundaries. 

• V is the space of elementwise- constant vectors with discontinuities across el- 
ement boundaries permitted. 

Example 2.9 (P1 DG -P2). The finite element pair known as P1 DG -P2 J2f is 
defined as follows: 

• The mesh M. is composed of triangular elements. 

• H is the space of elementwise- quadratic functions h s which are continuous 
across element boundaries. 

• V is the space of elementwise-linear vectors with discontinuities across ele- 
ment boundaries permitted. 

Each of these examples satisfy both conditions in Definition [23} condition 1 holds 
because taking the gradient of a elementwise polynomial n — 1 which is continuous 
across element boundaries results in a vector field which is discontinuous across ele- 
ment boundaries and is composed of elementwise polynomials of one degree n, and 
condition 2 holds since the velocity space uses the same continuity constraints for 
normal and tangential components e.g. both components are allowed to be discon- 
tinuous. This defines a whole sequence of high-order PnDG _ P( n +l) element pairs. 
Similar elements can be constructed on quadrilateral elements. Since we only require 
these two conditions to prove our optimal balance property which holds on arbitrary 
meshes, we can also construct finite clement spaces on mixed meshes composed of 
quadrilaterals and triangles, for example. It is also possible to use p-adaptivity in 
which different orders of polynomials are used in different elements, as long as the 
conditions are satisfied. 



3. Geostrophic balance properties. In this section we prove that when finite 
element pairs which are chosen to satisfy both of the conditions in Definition |2.5| 
their discrete geostrophically-balanced velocities satisfy the discrete divergence-free 
condition exactly: this means that the discrete geostrophically-balanced states are 
exactly steady states of equations ( 2.5||2.6 1. This is the result of Theorem 3.3 in this 
section, which makes use of conditions 1 and 2 of Definition 2.5 We also prove an 



inf-sup condition for these finite clement pairs (making use of condition 1 of Definition 
2.51 which provides a lower bound (independent of edge-lengths in the mesh) for the 
discrete gradient operator applied to non-constant functions in H . This lower bound 
prohibits the existence of spurious pressure modes which render a finite element pair 
unsuitable for geophysical flow problems; it also allows one to prove the convergence 
of the numerical solutions at the optimal rate obtained from approximation theory. 
This condition is stated and proved in Theorem |3.7| 



3.1. Optimal geostrophic balance. We first prove the following lemma which 
illustrates the embedding properties of our family of finite clement pairs. 

Lemma 3.1 (Pointwise gradient le mma) . Let (iZ, V) be a finite element pair 



chosen to satisfy condition 1 in Definition 2.5 Let Let q G V by the discrete gradient 
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oft] 6 G H defined by 

r w s ■ q 5 d V = I w s -VrfdV, (3.1) 



for all test functions w s G V . Then q s is the pointwise (strong) gradient of rf defined 
by 

q s (x) = Vr] S (x), Vcc G ft. 

Proof. Since condition 1 is satisfied, we may choose a test function 

w s = q s - Vrf G V, 



and substitution into equation 3.1 gives 



0= / w s ■ (q 5 - Vn s )dV 

Q 

S 



\q° - V77 a | dV 

Since the L2-norm only vanishes for elements of H if they are identically zero, we 
conclude that q s — \7-q 6 as required. □ This lemma appears at first sight to be a 
tautology but since the discrete gradient q s can be th ought of as the ^-projection of 



Vrj S into V, it requires condition 1 of Definition |2.5| to be satisfied, and it is not the 



case for equal-order element pairs such as Pl-Pl, for example. We shall make use of 



this lemma in proving Theorem 3.7 



The following lemma extends this technique to show that if the discrete geostrophic 
balance relation is satisfied by functions taken from a finite clement pair in our family 
of discretisations, then the exact geostrophic balance condition is actually satisfied at 
each point. 

Lemma 3.2 (Embedding lemm a). L et (H,V) be a finite element pair chosen to 



satisfy both conditions in Definition 2.5 Let u € V and r) < £ H satisfy the discrete 
geostrophic balance relation 

f [ w s -u 5± dV=-g [ w s -\7 v s dV, (3.2) 
Jn Jn 

,.s 



for all test functions w G V. Then 

= -gVr) 5 {x), Vx e SI. (3.3) 



Proof. Following the technique of the previous lemma, we note that conditions 1 
and 2 mean that we may choose a test function 

w s = fu s± + g\7r] s G V, 
and substitute into equation |3.2| to obtain 

= J w s ■ (fu 5± + g\/n s ^j d V 
fu s± + gX7r] 5 2 dV 

n 

l/^+sVr?- 5 !!?, 
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hence the result. □ 

This lemma is useful for initialising the system variables in a balanced state 
since the geostrophic balance relation can be evaluated pointwise instead of requiring 
a projection to be computed. Next we apply this lemma to prove the following 
theorem, which states that finite element pairs have the optimal property that these 
geostrophically balanced states are steady solutions of equations ( 2.5||2.6 ). 



Theo rem 3.3 (Steady geostrophic states). Let u° € V and rf = H satisfy 



equation {3.2), and let dfl be a contour for rf (so that the balanced velocity field 



obtained from equation 3.3 satisfies the boundary condition). Then u 
steady solutions of equation (2.5 2.6). 



and rf 



Proof. Substitution into equation (2.51, and choosing w = u gives 



dt 



u s \ 2 dV = 0, 



and hence uf — 0. It remains to show that 



rfdV = D Vf -uMV =0, 
Jn 

v ' 

divergence integral 



(3.4) 



for all test functions 6 s G H. 



By lemma 3.2 equation 3.3 is satisfied, which we may substitute into the diver- 
gence integral to obtain 



V-VdV. 



(3.5) 



The right-hand side integral in this equation can be shown to vanish for all (f> s and r] S 
in H 1 (which contains our velocity space H): the proof is obtained by taking 4> s and 
r] S as the limit of a convergent sequence of continuous functions in H 1 for which the 
sequence can be shown to vanish after integration by parts (see [8], for example). Here 
we provide a more direct proof which is obtained by integrating by parts separately 
in each element (since the gradients are discontinuous across element boundaries) to 
obtain 



V(j> 5 ■ VV d V 



{integration by parts} 



E 

EeM 



• V-VdV 



l5 V ■ vVdF 



<j) s n ■ V-VdV 




4 VVlldV 



TeM,rndn=% 



where E indicates an element in the mesh M. , dE is the boundary of E with outward 
pointing normal n, T indicates an edge in the mesh M. with an arbitrary chosen 
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normal direction, and [[u]] indicates the jump in the normal component of a vector 
field u across an edge T. In the final line, the first term vanishes since <fi is continuous 
across element boundaries, and the tangential components of V(j) S are also continuous 
(to see this, note that the tangential gradient of <j) 8 on the boundary is obtained by 
integrating in the direction of the boundary where r] 5 is continuously-differentiable) . 
The second term vanishes since the domain boundary 9f2 is a contour for rj S (from 
the assumption in the theorem), and so the tangential derivative of rf vanishes there. 



Hence, equation (|3.4| is satisfied. Choosing <p d — rf yields 



d 

dt 



6 

n 



r, 6 ) dV = 0, 



and hence r/f = and we have a steady state solution of equations ( 2.5||2.6 l. □ 
Remark 3.4. This proof corrects the proof presented in EV. 
Remark 3.5. Note that this theorem does not depend in any way on the structure 
of the mesh M. and so applies to arbitrary finite element discretisations on unstruc- 
tured meshes, provided that the conditions of Definition \2.5] are satisfied. 

We checked this theorem numerically by taking a completely unstructured mesh, 
randomly generating r] 5 and u s fields which satisfy the conditions of the theorem, 
and integrated equations ( 2.5j|2.6 1 using the implicit-midpoint rule. The problem 



was solved in dimensionless variables in a 1 x 1 square domain with Rossby number 
Ro = 0.1 and Froude number Fr = 1.0, with a timestep size At = 0.01 for 1000 steps. 
The maximum relative error between the initial and final fields was numerically 
zero {i.e. round-off error was observed) for each random realisation over hundreds of 
tests. Some example fields are shown in figure |3.1| These images illustrate that the 
optimal balance result is completely independent of the mesh and the smoothness of 
the solution. Some more general convergence tests using Kelvin waves are provided 
in g]. 

3.2. Inf-sup condition: pressure modes. We next prove that the discretisa- 
tion given by equations ( 2.5|2.6 1 using a finite element pair satisfying our two condi- 



tions does not have any spurious pressure modes. These are modes which converge (in 
the limit as the mesh edge-lengths go to zero) to steady solutions with zero velocity 
even though the free surface displacement if is not flat. These modes are catas- 
trophic for discretisations of the nonlinear equations since they can grow to dominate 
the solution when coupled to the physical modes. These modes are not present if it 
is possible to bound the gradient of non-constant functions in H from below as the 
edge-lengths go to zero. This bound is expressed in the following inf-sup condition: 

Definition 3.6 (Inf-sup condition). The inf-sup condition for a finite element 
pair (H, V) requires that 

sup Jn > J / {^Y&V, V/ € H, (3.6) 

"' eV yjla \ w5 \ 2dv V Ja 

where (3 is independent of the edge lengths in the mesh M. . 

Theorem 3.7 (Inf-sup theorem). Let (H,V) be a finite element pair which 



satisfies condition 1 in definition 2.5 Then [H,V) satisfies the inf-sup condition in 
definition\3.6\ 



Proof. The supremum in equation (3.6 1 is bounded from below by any particular 
choice of test function w s . By condition 1 X7<fi s e V and so we may choose w s = 
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J.142 0.284 0.426 0.566 0.710 0.853 0.995 



Fig. 3.1 

in Theorem 



3.3 



Figure showing random rf fields used to verify the optimal balance property described 
Randomly generated n s and u 6 fields ( constrained to satisfy the balance conditions 



of the theorem) were used as initi al condit ions for a numerical integration with the implicit-midpoint 
rule in time applied to eguations \2.5\2.6\ . The difference between the initial and final r) S fields was 
observed to be numerically zero for all random realisations that were tested. Note that the steady 
state does not depend on the smoothness of the solution or on the mesh structure. 



Substitution then gives 



sup /tX-V^d^ / n V^-V^dF 





\ 2 dV 


/ / W 


2 dV 


/ Ju 






2 dV 



b s2 dV, 



> VK^\ / 4> s2 dV, 
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where A m i n is the minimum non-zero eigenvalue of (minus) the Laplacian on f2, having 
made use of the Rayleigh quotient 



mm 



*7o j;,.v-\ir ' 
□ 

We note that the family of finite element pairs considered in this paper corre- 
sponds to one particular option for discrete differential forms satisfying a discrete 



de Rham complex condition described in [2J. In fact, condition 1 of Definition 2.5 
corresponds to all of the options described in that paper. 

A consequence of the inf-sup theorem, as discussed in O Q] , is that the solutions 
of the wave equation converge at the optimal rate described by the theory of numerical 
interpolation, as described in the following corollary: 

Corollary 3.8. Given an interval [0, T], there exists a constant C(T) such that 

\\u(;T)-u 5 (;T)\\ L2 + h(.,T) - V S (;T)\\ L2 < C{T)h k+l 



where (u,n) is the solution of equat ions ^2.1^2^2 ) with initial conditions (u ,r] ), 



(u s ,i] S ) is the solution of equations ( 2.^2^6 ) with initial conditions (mq,?7q) which 
satisfy the interpolation condition 

ll«o-«olU a + \\Vo-Vo\\l 2 < ch k+1 

and where k is the minimum of the orders of the elementwise polynomials used to 
construct u s and i] S . 

Proof. The proof follows using standard mixed finite element techniques, namely 
obtaining a bound on the L2-norm of the solution variables u s and if which requires 
the inf-sup condition. See [IT], [5] or [JJ, for example. □ 

For the P1dg _ P2 element pair, this convergence property (in this case 2nd- 
order convergence since k = 1) was confirmed from numerical experiments for the 
2-dimcnsional wave equation in [5], and for the rotating linear shallow-water equa- 
tions on an /-plane in [4]. 

4. Three dimensional incompressible flow. In this section we briefly dis- 
cuss the extension of these properties to the equations of three dimensional rotating 
stratified nonhydrostatic incompressible flow (Boussinesq equations) which, together 
with their hydrostatic counterparts, are used in ocean modelling. 

4.1. Model equations and geostrophic states. The full equations of motion 

are: 

u t + (u ■ V)u + fk x u = - — X7p+ kb, (4.1) 

Po 

V-w = 0, (4.2) 
T t + u ■ VT = 0, (4.3) 

where u is the three-dimensional velocity field, / is the Coriolis parameter, k is the 
unit upward vector, p is the pressure, b is the buoyancy and T is the temperature 
(salinity would also be included in a full ocean model but does not add anything to 
this discussion) . The system is closed by specifying an equation of state in which the 
buoyancy b is defined as a function of T and p. Here we have made the traditional 
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approximation which restricts the rotation vector to the vertical axis, and the /-plane 
approximation for which / is a constant. 

In three dimensions there are two geophysical balances which are present in slowly- 
varying large-scale flows which arise when the acceleration terms are small compared 
to the Coriolis and buoyancy. In the vertical direction we obtain the hydrostatic 
balance 

- — Pz + b = 0, (4.4) 
Po 

which can be imposed as a constraint in a hydrostatic model, allowing the pressure to 
be computed explicitly from the buoyancy by vertical integration. In the horizontal 
direction we again obtain the geostrophic balance 

11 . . 

u = fP y , v = — -p x . (4.5) 

Pof Pol 

It is simple to check again that V • u = Q for these balanced states. In the case of 
incompressible flow, this means that the balanced states are solutions of the equations 
given above. 

One can again take a numerical discretisation of equations ( 4.1|4.3 1, construct the 



discrete balanced solutions and check if the discrete form of equation |4.2| is satisfied 
exactly. If there is some residual, then this means that it is not possible for exactly- 
balanced states to exist in the numerical discretisation, which can lead to to the 
generation of spurious internal inertia-gravity waves. For example, if a pressure- 
projection method is used for timestepping (as is typical for non-hydrostatic models) 
then the time-integrator has two stages: the first stage takes a momentum step using 
the pressure from the previous timestep, then the solution is projected back to satisfy 
the discrete form of equation |4.2| If the balanced states do not satisfy this equation, 
then each projection will generate further spurious unbalanced motion. 

In this section we shall briefly describe an extension of our family of finite- clement 
pairs to the three-dimensional case, and describe the extension of our optimal balance 
results. Since we are still concerned with wave propagation we linearise equations 
( 4.1|4.3| about the state 



to obtain 



0, p z =p b, T = T(z), 



u t + fkxu = — -Vp' + kjT 1 , (4.6) 
Po 

V • u = 0, (4.7) 

T' t + u 3 T z = 0, (4.8) 



where 7 is a suitable positive constant. We shall drop the primes for the rest of the 
section. For these equations, steady balanced states given by 

"3 = 0, p z =p 'fT, u = -—p y , v=-^-p Xl 

Pof Pof 

satisfy V • u — and hence are admissible solutions of the equations. 
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4.2. Finite element formulation. Defining a Galcrkin finite clement method 
for these equations requires us to define a finite element space for the temperature 
variable T s . The choice of this space is not very important for the discussion of 
geostrophically-balanced states, so it will not be developed much here, except to note 
that it is important to ensure that there are sufficiently many states which satisfy 
a discrete hydrostatic balance, otherwise the representation of hydrostatic balance 
will be poor, leading to spurious non-hydrostatic motion. The velocity space V and 
the pressure space H are defined by the three dimensional extensions of definitions 

IH 

We now define the Galerkin finite element method for equations ( 4.6|4.8 1 as fol- 
lows: 

Definition 4.1 (Galerkin finite element method for 3D wave propagation). 
u s (x,t), p s (x 7 t), T s (x,t) are the semi- discrete solutions of the Galerkin finite el- 



ement discretisation of (4-6 4-8) if 



u s (;t)eV, P s (-,t)eH, T s (-,t)eQ, Vte[0,T], 



and 
d 

d7 



_d 
dt 



w d -u d dV + f / w d -kxu d dV = / w" 

Jn Pa Jn 



9 5 T 5 d V 



Vcf) d ■ u d dV 



fu s Tt d V 



0. 

0, 



Vp d dV + -fk- I inTplH) 
(4.10) 
(4.11) 



for all test functions w 6 V, (f> € H, 9 6 G. We again require our element pair 



(H, V) to satisfy the conditions in Definition 2.5 extended to three dimensions (with 



the _L operator replaced by the fcx operator). We next define the discrete geophysical 
balances using this discretisation. 

Definition 4.2 (Discrete hydrostatic and geostrophic balance). The solution 
variables u s , p s and T s satisfy the hydrostatic and geostrophic balances if 



kxu d dV 



1 



Pa Jn 



w d ■ Vp A d V + 7fc ■ 



w s T s d V. 



(4.12) 



The vertical component specifies hydrostatic balance and the horizontal component 
specifies geostrophic balance. This definition allows us to state the optimal bal- 
ance theorem for three-dimensional incompressible flow, which we give in the next 
subsection. 

4.3. Optimal balance properties, inf-sup theorem and optimal pressure 
matrix property. Theorem 4.3 (Optimal balance for three-dimensional incom- 
pressible flow). Let u s € V , p s £ H and T s € O be chosen so that equation {^.12) is 
satisfied, with zero vertical velocity 
condition 



and the pressure p satisfying the pointwise 



T H -Vp s = o, Vxedn, 



where tjj is the horizontal tangent vector to dfl so that the boundary is a streamline for 
the balanced flow (consistent with the boundary condition for u s ) . Then (u s ,p 5 ,T S ) 



is a steady state solution of equations (4.9 J^.ll) 
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Proof. It is simple to check that T s t — 0, uf — by inserting the solutions into 
the equations and noting that all of the terms vanish, and it remains to check that 
equation |4.10| is satisfied. The proof proceeds exactly as the proof of Theorem |3.3| by 



first noting that the two conditions in Definition 2.5 imply that 

J Pa 

pointwise, which we then insert into equation |4. 10| to obtain 

V0 s -u s dV= f </»^-^dy = 0, 
Jn 

using a similar argument to the previous proof. □ 

This means that balanced states are exactly steady and do not generate any 
spurious internal waves. Note that this theorem is again completely independent 
of the mesh structure and so the finite element pairs in this family are ideal for 
representing balanced flows on unstructured meshes such as those proposed in [15] , 

As for the 2D case, it is necessary for the finite element spaces to satisfy an inf- 
sup condition, so that solutions of the linearised equations converge at the optimal 
rate defined by approximation theory. In the case of incompressible flow, one also 
forms a pressure Poisson equation by composing the discrete divergence and gradi- 
ent operators, and if the inf-sup condition is not satisfied then there are very small 
spurious eigenvalues in the matrix which make iterative solvers very slow (see |10] , 
for example). Our family of finite clement pairs satisfy the inf-sup condition in three 
dimensions which can be shown by simple extension of the proof of Theorem |3.7| For 
incompressible flow, we shall use the same techniques to prove further properties of 
the discrete Poisson matrix. 

For the continuous equations, the Poisson equation is formed by taking the diver- 



gence of equation (4.6 1 and applying equation (4.7 1 to obtain 

V 2 p = Vr = -(p kjT- PQ fkxu), (4.13) 

which specifies an equation for p given T and u. This equation must be solved at each 
timestep to calculate the pressure field. When the Galerkin finite element method is 
applied, this specifies a coupled system of equations for p given by 



(4.15) 
(4.16) 



[ w s 


r 5 dV = 


Jn 




[ w s 


q s dV = 


Jn 




f 


q s dV = 


Jn 





V<f> 5 ■ r s < 
n 



for p s € H, q 5 , r 5 £ V and for all test functions <p s € H and w s 6 V. In practise, the 
variables q s and r s are eliminated to obtain an equation for p s . This then ensures 



that equation (4.10) is satisfied at each instance in time (or each timestep, having 
discrctised the equations in time). This system gives rise to a matrix equation for the 
basis function coefficients which, in general can have a larger sparsity pattern since 
it involve the product of several matrices (see |10j . for example). In the following 
theorem, we show that this sparsity pattern is reduced when condition 1 of Definition 
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|2.5| is satisfied, which is a further useful property of the family of finite clement pairs 
discussed in this paper. 

Theorem 4.4 (Optimal sparsity of pressure matrix). Let p s , q s , r s be the 
solutions of equations (4--l^\^.15 ). The 



len 



V(t>° ■ Vp° d V = - / V0" • r d V, (4.17) 
o Jn 



for all test functions (f> s € H , which is the usual Galerkin discretisation for equation 



(4-13) 



Proof. We again make use of lemma [3TT] which states that 

q S = Vp 5 



at each point. Substitution into equation (4.16) gives the result. □ Since the usual 



Galerkin discretisation of equation (4.131 results in a single equation that does not 
require the elimination of variables, this results in a much sparser stencil for the matrix 
equation for the basis function coefficients of p . 

5. Summary and outlook. In this paper we defined a large family of finite 
clement discretisations for the rotating shallow-water equations and the three dimen- 
sional equations of rotating stratified incompressible flow. When applied to the linear 
rotating shallow water equations, these discretisations were shown satisfy the optimal 
property that geostrophically-balanced states are completely steady, which mirrors 
a property of the solutions of the continuous equations. It was also shown that the 
discretisations in the family satisfy an inf-sup condition which prohibits the existence 
of spurious pressure modes. This makes the discretisations in the family strong can- 
didates for use in NWP and ocean modelling in cases where triangular elements are 
required, e.g. to allow the use of geodesic grids. Furthermore, the proofs are inde- 
pendent of the mesh structure which means that the family of discretisations produce 
stable results in the presence of adaptive mesh refinement. We then discussed the 
extension of the family to three-dimensional incompressible flow, required for ocean 
modelling, and showed that the discretisations in the family result in exactly steady 
balanced states on completely arbitrary unstructured meshes in three dimensions. In 
addition, the properties of the family were used to show that the matrix obtained 
from the discretised pressure Poisson equation has an optimally sparse stencil. 

In future work, we will develop and test discretisations of the fully-nonlinear 
equations using choices from our family of finite clement spaces, particularly from 
the PnDG _ P( n +l) sub- family applied to the rotating shallow- water equations and 
the three dimensional incompressible equations. Elements of that sub-family have 
discontinuous velocity, which allows a discontinuous Galerkin treatment of the ad- 
vection terms in the momentum equation. We also plan to compute numerical dis- 
persion relations for these discretisations when applied to the geodesic grid, in order 
to make comparisons with other element pairs and finite difference methods on this 
grid. The P1dg _ P2 element is currently being implemented in the Imperial College 
Ocean Model, and will be developed into a p-adaptive scheme in which different order 
polynomials are used in different elements. 
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